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ABSTRACT 

A coupled finite element /boundary element capability is described for 
calculating the sound pressure field scattered by an arbitrary submerged 3-D 
elastic structure. Structural and fluid impedances are calculated with no 
approximation other than discretization. The surface fluid pressures and 
normal velocities are first calculated by coupling a NASTRAN finite element 
model of the structure with a discretized form of the Helmholtz surface 
integral equation for the exterior fluid. Far-field pressures are then 
evaluated from the surface solution using the Helmholtz exterior integral 
equation. The overall approach is illustrated and validated using a known 
analytic solution for scattering from submerged spherical shells. 


INTRODUCTION 

Two fundamental problems of interest in structural acoustics are the 
calculation of the far-field acoustic pressure field radiated by a general 
submerged three-dimensional elastic structure subjected to internal time- 
harmonic loads and the calculation of the far-field acoustic pressure 
scattered by an elastic structure subjected to an incident time-harmonic wave 
train. These problems are usually solved by combining a finite element model 
of the structure with a fluid loading computed using either finite element 
[1-3] or boundary integral equation [4-11] techniques. 

Although both approaches are computationally expensive for large 
structural models, the fluid finite element approach is burdened with the 
additional complications caused by the approximate radiation boundary 
condition at the outer fluid boundary, the requirements on mesh size and 
extent, and the difficulty of generating the fluid mesh [1,3]. 
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In contrast, the boundary integral equation (BIE) approach for generating 
the fluid loading is mathematically exact (except for surface discretization 
error) and requires little or no additional modeling effort to convert an 
existing model of a dry structure for use in submerged analyses. The saving 
in engineering time, however, is partially offset by the somewhat greater 
computing costs associated with the BIE approach. 

Several general BIE acoustic radiation capabilities have been developed 
previously [4, 7, 11]. One, called NASHUA [11], couples a NASTRAN finite 
element model of a dry structure with a fluid loading calculated by a 
discretized form of the Helmholtz surface integral equation. NASHUA is the 
only capability developed for a widely-used, public domain, general purpose 
structural analysis code. Here we present an extension to NASHUA to handle 
also the problem of acoustic scattering from general three-dimensional elastic 
structures. 

The primary purposes of this paper are to summarize the theoretical 
basis for NASHUA and to demonstrate its validity for scattering by showing 
results of calculations for the elementary problem of plane-wave scattering 
from a thin spherical shell. 


THEORETICAL APPROACH 

Consider an arbitrary submerged 3-D elastic structure subjected to either 
internal time-harmonic loads or an external time-harmonic incident pressure 
wave train. The matrix equation of motion for the structural degrees of 
freedom (DOF) can be written as 


Zv = F - GAp (1) 


where Z = structural impedance matrix (dimension s x s), 

v = complex amplitude of the velocity vector for all structural DOF 
(wet or dry) in terms of the coordinate systems selected by the 
user (s x r), 

F = complex amplitude of the vector of mechanical forces applied to the 
structure (s x r), 

G = rectangular transformation matrix of direction cosines to transform 
a vector of outward normal forces at the wet points to a vector of 
forces at all points in the coordinate systems selected by the user 
(s x f), 

A = diagonal area matrix for the wet surface (f x f), and 
p = complex amplitude of total pressures (incident + scattered) applied 
at the wet grid points (f x r). 

In this equation, the time dependence exp(ia)t) has been suppressed. In the 
above dimensions, s denotes the total number of structural DOF (wet or dry), 
f denotes the number of fluid DOF (the number of wet points), and r denotes 
the number of load cases. In general, surface areas in NASHUA are obtained 
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from the NASTRAN calculation of the load vector resulting from an outwardly 
directed static unit pressure load on the structure's wet surface. 

In Equation (1), the structural impedance matrix Z (the ratio of force 
to velocity) is given by 


Z = (-w 2 M + iwB + IC)/iu) 


( 2 ) 


where M, B, and K are the structural mass, viscous damping, and stiffness 
matrices, respectively, and w is the circular frequency of excitation. For 
structures with material damping or a nonzero loss factor, K. is complex. 

The total fluid pressure p satisfies the reduced wave equation 


V^p + k^p - 0 (3) 


where k * w/c is the acoustic wave number, and c is the speed of sound in the 
fluid. Equivalently, p is the solution of the Helmholtz integral equation 
[7,12] 


( p(x_* )/2 “ Px, x' on S 

f p(x)(3D(r)/3n)dS - / q(x)D(r)dS = { (4) 

s s ( p(x' ) - Pj, x' in E 


where S and E denote surface and exterior fluid points, respectively, px is 
the incident free-field pressure, r is the distance from x. to i’ (Figure 1), 


FLUID 


Pi 



Figure 1 - Notation for Helmholtz Integral Equations 
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D is the Green’s function 


D(r) = 


and 


q = 9p/3n = -ioipv n 


(5) 

( 6 ) 


where p is the density of the fluid, and v n is the outward normal component of 
velocity on S. As shown in Figure 1 , jc in Equation (4) is the position vector 
for a typical point Pj on the surface S, x' is the position vector for the 
point P-£ which may be either on the surface or in the exterior field E, the 
vector r^ - jc’ - jc, and ri is the unit outward normal at Pj » We denote the 
lengths of the vectors jc, x ’ , and jr by x, x f , and r, respectively. The normal 
derivative of the Green’s function D appearing in Equation (4) can be 
evaluated as 


3D(r)/3n 35 (e~ ilcr /4-rrr) (ik + 1/r) cos 0 


(7) 


where 0 is defined as the angle between the normal _n and the vector £, as 
shown in Figure 1 . 

The substitution of Equations (6) and (7) into the surface equation (4) 
yields 

p(x’)/2 - / p(x) (e”i^ r /47rr) (ik + 1/r) cos 0 dS 
S 

= imp / v n (x) (e“ ikr /4Trr )dS + p x (8) 

S 


where x’ is on S. This equation is an integral equation relating the total 
pressure p and normal velocity v n on S. If the integrals in Equation (8) are 
discretized for numerical computation (the details of which were presented 
previously [11]), we obtain the matrix equation 


Ep = Cv n + p x 


(9) 


on S, where p is the vector of complex amplitudes of the total pressure on the 
structure’s surface, E and C are fully-populated, complex, non-symmetric, 
frequency-dependent matrices, and pj is the complex amplitude of the incident 
pressure vector (if any). The number of unknowns in this system is f, the 
number of wet points on the fluid-structure interface. 
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The normal velocities v n in Equation (9) are related to the total 
velocities v by the same transformation matrix G: 


v n = ° Tv 


( 10 ) 


where T denotes the matrix transpose* If velocities v and v n are eliminated 
from Equations (1), (9), and (10), the resulting equation for the coupled 
fluid-structure system is 


Hp = Q + P! (11) 

where 

H = E + CGT Z _ 1 GA (12) 

and 

Q = CG t Z -1 F (13) 


Since H and Q depend on geometry, material properties, and frequency, Equation 
(11) may be solved to yield the total surface pressures p. The vector v of 
velocities at all structural DOF may be recovered by solving Equation (1) 
for v: 


v = Z _1 F - Z“1 GAp 


(14) 


Surface normal velocities v n may then be recovered by substituting this 
solution for v into Equation (10). The triangular factors of Z are saved 
when first generated in Equation (12) and reused in Equations (13) and (14). 

The free-field incident pressure for planar or spherical waves is 
calculated in the following way. Consider an arbitrary surface S subjected 
to a time-harmonic incident wave train of speed c and circular frequency cj as 
shown in Figure 2. The complex amplitude p-j- of the incident pressure at a 
typical point i is given by 


PI 


' p Q exp(i0^), plane wave 

p 0 expd^) s/(s - <5^), spherical wave 


(15) 


where p Q is the pressure amplitude at the coordinate origin (where the phase 
angle is arbitrarily set to zero), s is the length of the vector s_ (the vector 
from the origin to the source), 6 ^ is the distance (positive or negative) 
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from the origin to the wavefront through point i, and 9^ is the phase angle 
for point i. These last two quantities can be computed using 


0i = k 6 ± 


( 16 ) 


and 


I x ± * s/s, plane wave 

s - |xi - sj , spherical wave 


(17) 


where 2LL is the position vector for point i. For plane waves, only the 
direction cosines of are used. Positive values of 9x and 6-^ for a point i 
correspond to that point’s being closer to the source than the origin is. 

To summarize, the NASHUA solution procedure uses NASTRAN to generate K, 
M, B, and F and to generate sufficient geometry information so that E, C, G, 
A, and px can be computed by a separate program (SURF). Then, given all 
matrices on the right-hand sides of Equations (12) and (13), NASTRAN DMAP is 
used to compute H and Q. Equation (11) is then solved for the surface 
pressures p using a new block solver (OCSOLVE) written especially for this 
problem by E.A. Schroeder of the David Taylor Naval Ship R&D Center. Next, 
NASTRAN DMAP is used to recover the surface normal velocities v n and the 
vector v of velocities at all structural DOF. This completes the surface 
solution. In general, this approach combines in a highly automated fashion a 
finite element model of the structure with a Helmholtz boundary integral 
equation model of the fluid. 


/ 



SOURCE 


Figure 2 - Geometry for Calculating Incident Pressure 
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The Far-Field Solution 


Given the solution for the total pressures and velocities on the surface, 
the exterior Helmholtz integral equation. Equation (4), can be integrated to 
obtain the radiated (or scattered) pressure at any desired location x * in the 
field. We first substitute Equations (6) and (7) into Equation (4) to obtain 
a form suitable for numerical integration: 


p(x') = / s [iu)pv n (x.) + + 1 / r )p(2L^ cos 3] (e ilcr /4Trr)dS (18) 


where all symbols have the same definitions as were used previously, and x* 
is in the exterior field. Thus, given the total pressure p and normal 
velocity v n on the surface S, the radiated or scattered pressure at x 1 can be 
determined by numerical quadrature using Equation (18). 

In applications, however, the field pressures generally of interest are 
in the far-field, so we use instead an asymptotic (far-field) form [11] 
of Equation (18): 


p(x') = (ike“ ikx /47TX 1 ) / [pcv (x) + p(x) cos 3je ikx cos a dS (19) 

where a is the angle between the vectors x and x* (Figure 1), and, for far- 
field points, cos B is computed using 


cos 3 = n-x'/x' 


( 20 ) 


We note that, since Equation (19) is a far-field formula, the pressure varies 
inversely with distance x* everywhere so that any range x T may be used in its 
evaluation (e.g., 36 inches). 


OVERVIEW OF NASHUA SOLUTION PROCEDURE 

The overall organization and setup of the solution procedure is 
summarized in Figure 3. NASTRAN appears four times in the procedure; to 
distinguish one NASTRAN execution from another, the integers 1-4 are appended 
to "NASTRAN" in the figure. 

A separate NASTRAN model is prepared and run (Step 1 in Figure 3) for 
each unique set of symmetry constraints. Since up to three planes of 
reflective symmetry are allowed, there would be one, two, four, or eight such 
runs. Step 1 generates files containing geometry information and a checkpoint 
file for subsequent use in the other steps. 
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For each symmetry case and drive frequency, the Step 2 sequence is run 
in a single job. The SURF program reads the geometry file generated by 
NASTRAN in Step 1 and, using the Helmholtz surface integral equation, 
generates the fluid matrices E and C for the exterior fluid, the area matrix 
A, the structure-fluid transformation matrix G, the free-field incident 
pressure vector pj, and a geometry file to be used later by FAROUT (Step 3) 
for the field calculation. SURF is followed by a NASTRAN job which takes the 
matrices K, M, B, and F from Step 1 and the matrices E, C, A, and G from SURF 
and calculates H and Q according to Equations (12) and (13). Equation (11) 
is then solved for the surface pressure vector p by program OCSOLVE written 
by E.A. Schroeder of DTNSRDC. OCSOLVE is a general block solver for full, 
complex, nonsymmetric systems of linear, algebraic equations. The program 
was designed to be particularly effective on such systems and executes about 
20 times faster than NASTRAN’s equation solver, which was not designed for 
efficient solution of such systems of equations. NASTRAN is then re-entered 
in Step 2 with p so that the velocities v and v n can be recovered using DMAP 
operations according to Equations (14) and (10), respectively. 



NOTE: Each solid block is a separate job submission. 


Figure 3 — Summary of NASHUA Solution Procedure 
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The surface pressures, normal velocities, and full g-set displacements are 
then reformatted and merged into a single file (for each symmetry case) using 
program MERGE. 

Steps 1 and 2 are repeated for each symmetry case. After all symmetry 
cases have been completed and merged, program FAROUT (Step 3) is run to 
combine the symmetry cases and to integrate over the surface. FAROUT uses as 
input the geometry file generated by SURF (Step 2) and the surface solutions 
from the one, two, four, or eight symmetry cases generated by MERGE. The 
far-field pressure solution is obtained by integrating the surface pressures 
and velocities using the far-field form of the exterior Helmholtz integral 
equation, Equation (19). Output from FAROUT consists of both tables and 
files suitable for various types of plotting. 

The remaining steps in the NASHUA procedure are for graphical display. 
Deformed structural plots of the frequency response may be obtained by 
restarting NASTRAN (Step 4) with the checkpoint file from Step 1 and a results 
file from FAROUT. In addition, animated plots can be generated on the Evans 
& Sutherland PS-330 graphics terminal using the CANDI program (Step 6) written 
for the DEC/VAX computer by R.R. Lipraan of DTNSRDC [13]. If the rest of 
NASHUA is run on a computer other than the VAX, the NASTRAN UT1 file which is 
passed to CANDI must first be formatted (Step 4) for transfer to the VAX 
computer and then unformatted (Step 5) for reading by CANDI. 

X-Y plots of various quantities (both surface and far-field) versus 
frequency may be obtained using the general purpose interactive plotting 
program IPLOT [14] (Step 7). Polar plots of the far-field sound pressure 
levels in each of the three principal coordinate planes can also be generated 
using the interactive graphics program FAFPLOT [13] (Step 8) written by R.R. 
Lipman of DTNSRDC. 


FREQUENCY LIMITATIONS 

It is known that the fluid matrices E and C in the surface Helmholtz 
integral equation formulation are singular at the frequencies of the 
resonances of the corresponding interior acoustic cavity with Dirichlet (zero 
pressure) boundary conditions [3]. Although the NASHUA formulation described 
in the previous section was designed to avoid having to invert either E or C 
in Equations (12) and (13), the coefficient matrix H is also poorly conditioned 
at these frequencies (referred to as the "critical” or "forbidden” frequencies 
of the problem) [11]. Therefore, to be safe, the user should generally avoid 
excitation frequencies which exceed the lowest critical frequency for the 
geometry in question. 

For spheres, for example, the lowest critical frequency occurs at ka = it, 
where k is the acoustic wave number, and a is the radius. For long cylinders 
with flat ends, the lowest critical frequency occurs at ka ^ 2.4, where a is 
the radius. For short cylinders with flat ends, the lowest critical frequency 
is slightly higher than for long cylinders. 
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RESTRICTIONS ON MODEL 


Although the NASHUA solution procedure was designed to be general enough so 
that arbitrary three-dimensional structures could be analyzed, a few restrictions 
remain. In our view, however, none is a burden, since a NASTRAN deck for a dry 
structure modeled with low— order finite elements can usually be adapted for 
use with NASHUA in a few hours. The following general restrictions apply: 

1. All translational degrees of freedom (DOF) for wet points must be in 
NASTRAN’s "analysis set" (a-set), since all symmetry cases must have the same 
wet DOF, and the fluid matrices E and C involve all wet points. This 
restriction also affects constraints. Thus, constraints on translational DOF 
of wet points may not be imposed with single point constraint (SPC) cards, 
but must instead be imposed using large springs connected between the DOF to 
be constrained and ground. Generally, this restriction affects only those 
DOF which are constrained due to symmetry conditions. 

2. The wet face of each finite element in contact with the exterior 
fluid must be defined by either three or four grid points, since the numerical 
discretization of the Helmholtz surface integral equation assumes the use of 
low order elements. In particular, NASTRAN elements with midside nodes 
(e.g., TRIM6, IS2D8, or IHEX2) may not be in contact with the exterior fluid. 

3. Symmetry planes must be coordinate planes of the basic Cartesian 
coordinate system. 

4. No scalar points or extra points are allowed, since program SURF 
assumes that each point is a grid point. 

3. For cylindrical shells, the axis of the cylinder should coincide 
with one of the three basic Cartesian axes; for spherical shells, the center 
of the sphere should coincide with the basic origin. These restrictions 
facilitate the treatment of symmetry planes and the calculation of curvatures 
in program SURF. 

6. At least one degree of freedom in the model should be constrained 
with an SPC, MPC, or OMIT card so that the NASTRAN data block PL is generated. 

7. Thin structures with fluid on both sides should be avoided, since 
the formulations for the fluid matrices are singular if two wet points are 
coincident. A precise restriction is not known. 


TIME ESTIMATION 

On CDC computers, most of the computer time required to execute the 
entire NASHUA procedure is associated with the forward/backward solve (FBS) 
operation in Step 2, Equation (12), in which the matrix Z“*GA is computed 
given the triangular factors of Z and the matrix GA. Z is a complex, 
symmetric, banded matrix of dimension s x s, where s is the number of 
structural DOF in the problem, and GA is a real, sparsely-populated. 
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rectangular matrix of dimension s x f, where f is the number of fluid DOF 
(the number of wet points on the surface). This FBS time is proportional to 
the product of s , f, and W avg (the average wavefront for the stiffness matrix 
K), and, for large jobs, accounts for a substantial part (perhaps two-thirds) 
of the total time to make a single pass through the NASHUA Step 2 procedure. 

For example, consider a problem with the following characteristics: 

s = 2973 (number of structural DOF) 
f = 496 (number of fluid DOF) 

W av g = 129 (average wavefront of stiffness matrix) 

On the CDC Cyber 176 computer at DTNSRDC, the computer time ( ’’wall-clock" 
time) required to solve this problem in a dedicated computer environment for 
a single symmetry case and one drive frequency was about 30 minutes, of which 
19 minutes were spent in the FBS operation. 


PROGRAM VALIDATION 

For radiation problems, NASHUA has been validated previously [11] for 
submerged spherical shells driven internally by both uniform and non-uniform 
(sector) pressure loads. Here we demonstrate NASHUA'S ability to solve 
scattering problems by solving the problem of the submerged thin spherical 
shell subjected to an incident time-harmonic planar wave train, as shown in 
Figure 4. The solution of this problem exhibits rotational symmetry about 
the spherical axis parallel to the direction of wave propagation. The 
benchmark solution to which the NASHUA results will be compared is a series 
solution published in the Junger and Feit book [16]. 


STEEL 



WATER 


PLANAR 

WAVETRAIN 


Figure 4 - Plane Wave Scattering from an Elastic Spherical Shell 
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The problem to be solved has the following characteristics [17]: 


a ** 5 m 
h = 0.15 m 
E = 2.07 x 10H Pa 
v = 0.3 

p = 7669 kg/m^ 

T1 = 0 

p = 1000 kg/m^ 

c - 1524 m/s 


(shell radius) 

(shell thickness) 
(Young's modulus) 
(Poisson's ratio) 
(shell density) 

(shell loss factor) 
(fluid density) 

(fluid speed of sound) 


One octant of the shell was modeled with NASTRAN's CTRIA2 membrane/bending 
elements as shown in Figure 5. With 20 elements along each edge of the 
domain, the model has 231 wet points and 1263 structural DOF. Since the 
incident loading does not exhibit three planes of symmetry, the NASHUA 
solution of this problem requires decomposing the solution into both symmetric 
and antisymmetric parts of the problem, thus providing a good check on 
NASHUA'S ability to combine symmetry cases for scattering problems. 

The NASHUA model was run for 15 different drive frequencies in the 
nondimensional frequency range ka = 0.5 to ka = 5.0, where a is the shell 
radius. Two of the excitation frequencies are near a critical frequency. 

(The first 13 critical frequencies are located at ka = tt , 4.49, 5.76, 2tt , 

6.99, 7.73, 8.18, 9.10, 9.36, 3tt , 10.4, 10.5, and 10.9 [18].) Figure 6 shows 
a comparison between the NASHUA calculations and the series solution for the 
far-field scattered pressure in the forward direction (9 = 180 degrees). 



Figure 5 - Finite Element Model of One Octant of Spherical Shell 
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The ordinate of this figure is the normalized pressure |p r r/p 0 a|, where p r is 
the far-field scattered pressure at distance r from the origin, and p Q is the 
magnitude of the incident pressure. Clearly, the NASHUA solution agrees very 
well with the exact (i.e., converged series) solution. We note that the 
NASHUA calculation at ka 888 4.5 is adversely affected by the forbidden 
frequency at ka * 4.49. (Often, the effect of a forbidden frequency is very 
severe [11].) For several of the excitation frequencies, we also tabulate on 
the next page the far-field scattered pressure patterns. Again the agreement 
between the NASHUA calculations and the series solution is excellent, even at 
ka * 1.6, which is near a resonant peak. At the sharper resonant peaks, the 
results would be much more sensitive to small changes in frequency. 


DISCUSSION 

A very general capability has been described for predicting the acoustic 
sound pressure field scattered by arbitrary three-dimensional elastic 
structures subjected to time-harmonic incident loads. Sufficient automation 
is provided so that, for many structures of practical interest, an existing 



Figure 6 - Normalized Far-Field Pressure |p r r/p 0 a| Scattered in the Forward 
Direction (9 = 180 degrees) by a Spherical Shell; Solid Curve is Converged 
Series Solution [16], and Square Boxes are NASHUA Solution. 
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Table - Comparison of NASHUA Solution with Converged Series 
Solution for Scattering from Spherical Shell 


ka 

Angle 9 
(degrees ) 

Normalized Far-Field Pressure, Tpr r /Pn a l ! 

NASHUA 

Exact 

% Error 


0* 

0.0083 

0.0081 

2.5 


30 

0.0144 

0.0143 

0.7 


60 

0.0299 

0.0299 

0.0 

0.5 

90 

0.0481 

0.0481 

0.0 


120 

0.0626 

0.0626 

0.0 


150 

0.0708 

0.0708 

0.0 


180 

0.0734 

0.0733 

0.1 


0* 

0.0892 

0.0903 

1.2 


30 

0.0382 

0.0389 

1.8 


60 

0.0886 

0.0886 

0.0 

1.0 

90 

0.1926 

0.1930 

0.2 


120 

0.2208 

0.2210 

0.1 


150 

0.1893 

0.1887 

0.3 


180 

0.1662 

0.1652 

0.6 


0* 

3.146 

3.149 

0.1 


30 

1.993 

1.995 

0.1 


60 

0.325 

0.320 

1.6 

1.6 | 

90 

1.515 

1.498 

1.1 


120 

0.561 

0.540 

3.9 


150 

2.058 

2.092 

1.6 


180 

3.213 

3.245 

1.0 


0* 

0.166 

0.160 

3.8 


30 

0.105 

0.101 

4.0 


60 

0.068 

0.069 

1.4 

4.0 

90 

0.274 

0.269 

1.9 


120 

0.562 

0.554 

1.4 


150 

1.769 

1.757 

0.7 


180 

3.228 

3.205 

0.7 


* 9*0 corresponds to the back-scattered direction. 


NASTRAN structural model can be adapted for NASHUA acoustic analysis within a 
few hours. 

One of the major benefits of having NASHUA linked with NASTRAN is the 
ability to integrate the acoustic analysis of a structure with other dynamic 
analyses. Thus the same finite element model can be used for modal analysis, 
frequency response analysis, linear shock analysis, and underwater acoustic 
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analysis. In addition, many of the pre- and postprocessors developed for use 
with NASTRAN become available for NASHUA as well. 

The principal area in which NASHUA could be improved would be to remove 
the frequency limitation caused by the presence of the critical frequencies 
inherent in the Helmholtz integral equation formulation. With the limitation 
cylindrical shells, for example, can be safely analyzed by NASHUA only for 
ka < 2.4, where a is the radius. Since for some problems, it would be of 
interest to treat higher frequencies, the limitation should be removed. A 
conversion to a different formulation (e.g., Burton and Miller [8] or Mathews 
[10]) is being considered. 
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